%% Ergasia1_splines_parembolh_12.03.2009
%% Onoma: Georgiou Dimitrios, A.E.M.:569

%% Ergasia 1a
clc
clear
format compact
disp('Kwdikas gia mia kampylh fysikhs kyvikhs paremvolhs');
%% (1) Orismow twn dedomenwn gia tis meres ths evdomadas kai gia tous 2 deiktes
x = [1 2 3 4 5];    % meres th evdomadas
y = [7607.40 7544.60 7561.20 7526.40 7515.30];   % deikths gia TSE300
z = [1848.20 1836.50 1860.20 1843.40 1805.20];   % deikths gia NASDAQ

h = 1;  % h = x2 -x1 = x3 -x2 = … = xn+1 -xn

n = length(y)-1;  % n = vathmos tou polywnymou

%% (2) Arxikopoihs tou pinaka A gia opoiodhpote n kai prosdiorismos tou b

% Arxikopoihs tou pinaka A
A = 4 * eye(n-1,n-1);
for  i = 1:n-1
     if(i == 1) 
         A(i,i+1)= 1;
     elseif(i == n-1)
         A(i,i-1) = 1;
     else
         A(i,i-1)= 1;
         A(i,i+1)= 1; 
     end
end
A
% prosdiorismos tou dianysmatos b
% Symfwna me ton orismo dinw times sta dianysmata
for i=1:n-1
   b_y(i) =  (6/(h^2))* (y(i) - 2*y(i+1) + y(i+2)); % Dimiourgo ena dianysma by, n-1 orwn gia ton deikth TSE300
   b_z(i) =  (6/(h^2))* (z(i) - 2*z(i+1) + z(i+2)); % Dimiourgo ena dianysma bz, n-1 orwn gia ton deikth NASDAQ
end

%% (3) Ypologismos twn timwn mj, gia kathe synolo dedomenwn
% Sthn periptwsh ayth tha prepi na lisw to systhma Am=b

% Ypologizw ta m(2), m(3), ..., m(n) gia tou diktes TSE300 kai NASDAQ

tempy = A\b_y';  %deikth TSE300
tempz = A\b_z';  %deikth NASDAQ

% Ypologismos twn m(j), j=1,2,...,n,n+1. gia ton deikth TSE300
% Ypologismos twn m(j), j=1,2,...,n,n+1. gia ton deikth NASDAQ
for i=1:n+1
    if(i == 1)||(i == n+1)
        my(i) = 0;  %deikths TSE300
        mz(i) = 0;  %deikths NASDAQ
    else
        my(i)=tempy(i-1);   %deikths TSE300
        mz(i)=tempz(i-1);   %deikths NASDAQ
    end
end
my
mz

%% (4) ypologismos twn paramertwn ( a(j), b(j), c(j), d(j) ), twn
%% polywnymown kivikhs paremvolhs y = Sj(x)

ay = y(1:n);
az = z(1:n);
   
by = (y(2:n+1) - y(1:n))./h - h.*(my(2:n+1)+2*my(1:n))/6;
bz = (z(2:n+1) - z(1:n))./h - h.*(mz(2:n+1)+2*mz(1:n))/6;

   
cy = 0.5*my(1:n);
cz = 0.5*mz(1:n);

dy = (my(2:n+1) - my(1:n))./(6*h);
dz = (mz(2:n+1) - mz(1:n))./(6*h);

%% (5) Sxediasmos twn splines paremvolis sto diasthma[0,5] gia ta 2 synola
%% (6) kai Sxediasmos twn simiwn twn antoistixwn dedomenwn

xInt = 1 : 0.01 : 5;
for j = 1 : length(xInt)
	if xInt(j) ~= x(n+1)
        iInt(j) = sum(x <= xInt(j));
	else
		iInt(j) = n;
	end
end



yInt = ay(iInt) + by(iInt).*(xInt-x(iInt)) + cy(iInt).*(xInt - x(iInt)).^2 + dy(iInt).*(xInt - x(iInt)).^3;
plot(xInt,yInt,'g',x,y,'b*');   axis([0,5,7500,7950]);

pause;
zInt = az(iInt) + bz(iInt).*(xInt-x(iInt)) + cz(iInt).*(xInt - x(iInt)).^2 + dz(iInt).*(xInt - x(iInt)).^3;
plot(xInt,zInt,'g',x,z,'b*');   axis([0,5,1800,1870]);
% Oi times sta mesa ths tetarths
yInt = ay(iInt) + by(iInt).*(3.5-x(iInt)) + cy(iInt).*(3.5 - x(iInt)).^2 + dy(iInt).*(3.5 - x(iInt)).^3;
zInt = az(iInt) + bz(iInt).*(3.5-x(iInt)) + cz(iInt).*(3.5 - x(iInt)).^2 + dz(iInt).*(3.5 - x(iInt)).^3;
pause;

fprintf('I timh sta mesa ths tetarths:\n');
fprintf('Gia ton deikth TSE300 :');disp(sum(yInt)/length(yInt))
fprintf('\nGia ton deikth NASDAQ :');disp(sum(zInt)/length(zInt))

pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Ergasia 1b
clc
clear
format compact
disp('\nKwdikas me gia mia kivikhs spline paremvolhs me paravolikh simperifora sta akra');
%% (1) Orismow twn dedomenwn gia tis meres ths evdomadas kai gia tous 2 deiktes
x = [1 2 3 4 5];    % meres th evdomadas
y = [7607.40 7544.60 7561.20 7526.40 7515.30];   % deikths gia TSE300
z = [1848.20 1836.50 1860.20 1843.40 1805.20];   % deikths gia NASDAQ

h = 1;  % h = x2 -x1 = x3 -x2 = … = xn+1 -xn
n = length(y)-1;  % n = vathmos tou polywnymou

%% (2) Arxikopoihs tou pinaka A gia opoiodhpote n kai prosdiorismos tou b

% Arxikopoihs tou pinaka A
A = 4 * eye(n-1,n-1);
for  i = 1:n-1
     if(i == 1) 
         A(i,i) = 5;
         A(i,i+1)= 1;
     elseif(i == n-1)
         A(i,i) = 5;
         A(i,i-1) = 1;
     else
         A(i,i-1)= 1;
         A(i,i+1)= 1; 
     end
end
A
% prosdiorismos tou dianysmatos b
% Symfwna me ton orismo dinw times sta dianysmata
for i=1:n-1
   b_y(i) =  (6/(h^2))* (y(i) - 2*y(i+1) + y(i+2)); % Dimiourgo ena dianysma by, n-1 orwn gia ton deikth TSE300
   b_z(i) =  (6/(h^2))* (z(i) - 2*z(i+1) + z(i+2)); % Dimiourgo ena dianysma bz, n-1 orwn gia ton deikth NASDAQ
end

%% (3) Ypologismos twn timwn mj, gia kathe synolo dedomenwn
% Sthn periptwsh ayth tha prepi na lisw to systhma Am=b

% Ypologizw ta m(2), m(3), ..., m(n) gia tou diktes TSE300 kai NASDAQ

tempy = A\b_y';  %deikth TSE300
tempz = A\b_z';  %deikth NASDAQ

% Ypologismos twn m(j), j=1,2,...,n,n+1. gia ton deikth TSE300
% Ypologismos twn m(j), j=1,2,...,n,n+1. gia ton deikth NASDAQ
for i=2:n
        my(i)=tempy(i-1);   %deikths TSE300
        mz(i)=tempz(i-1);   %deikths NASDAQ
end
my(1) = my(2);  %deikths TSE300
mz(1) = mz(2);  %deikths NASDAQ

my(n+1) = my(n);  %deikths TSE300
mz(n+1) = mz(n);  %deikths NASDAQ
        
my
mz
%% (4) ypologismos twn paramertwn ( a(j), b(j), c(j), d(j) ), twn
%% polywnymown kivikhs paremvolhs y = Sj(x)

ay = y(1:n);
az = z(1:n);
   
by = (y(2:n+1) - y(1:n))./h - h.*(my(2:n+1)+2*my(1:n))/6;
bz = (z(2:n+1) - z(1:n))./h - h.*(mz(2:n+1)+2*mz(1:n))/6;

   
cy = 0.5*my(1:n);
cz = 0.5*mz(1:n);

dy = (my(2:n+1) - my(1:n))./(6*h);
dz = (mz(2:n+1) - mz(1:n))./(6*h);

%% (5) Sxediasmos twn splines paremvolis sto diasthma[0,5] gia ta 2 synola
%% (6) kai Sxediasmos twn simiwn twn antoistixwn dedomenwn

xInt = 1 : 0.01 : 5;
for j = 1 : length(xInt)
	if xInt(j) ~= x(n+1)
        iInt(j) = sum(x <= xInt(j));
	else
		iInt(j) = n;
	end
end
yInt = ay(iInt) + by(iInt).*(xInt-x(iInt)) + cy(iInt).*(xInt - x(iInt)).^2 + dy(iInt).*(xInt - x(iInt)).^3;
plot(xInt,yInt,'g',x,y,'b*');   axis([0,5,7500,7950]);

pause;
zInt = az(iInt) + bz(iInt).*(xInt-x(iInt)) + cz(iInt).*(xInt - x(iInt)).^2 + dz(iInt).*(xInt - x(iInt)).^3;
plot(xInt,zInt,'g',x,z,'b*');   axis([0,5,1800,1870]);

% Oi times sta mesa ths tetarths
yInt = ay(iInt) + by(iInt).*(3.5-x(iInt)) + cy(iInt).*(3.5 - x(iInt)).^2 + dy(iInt).*(3.5 - x(iInt)).^3;
zInt = az(iInt) + bz(iInt).*(3.5-x(iInt)) + cz(iInt).*(3.5 - x(iInt)).^2 + dz(iInt).*(3.5 - x(iInt)).^3;
pause;

fprintf('I timh sta mesa ths tetarths:\n');
fprintf('Gia ton deikth TSE300 :');disp(sum(yInt)/length(yInt))
fprintf('\nGia ton deikth NASDAQ :');disp(sum(zInt)/length(zInt))
pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Ergasia 1c
clc
clear
format compact
disp('Kwdikas me gia mia kivikhs spline paremvolhs me eniaia kivikh simperifora sta akra');
%% (1) Orismow twn dedomenwn gia tis meres ths evdomadas kai gia tous 2 deiktes
x = [1 2 3 4 5];    % meres th evdomadas
y = [7607.40 7544.60 7561.20 7526.40 7515.30];   % deikths gia TSE300
z = [1848.20 1836.50 1860.20 1843.40 1805.20];   % deikths gia NASDAQ

h = 1;  % h = x2 -x1 = x3 -x2 = … = xn+1 -xn
n = length(y)-1;  % n = vathmos tou polywnymou

%% (2) Arxikopoihs tou pinaka A gia opoiodhpote n kai prosdiorismos tou b

% Arxikopoihs tou pinaka A
A = 4 * eye(n-1,n-1);
for  i = 1:n-1
     if(i == 1) 
         A(i,i) = 6;
         A(i,i+1)= 1;
     elseif(i == n-1)
         A(i,i) = 6;
         A(i,i-1) = 1;
     else
         A(i,i-1)= 1;
         A(i,i+1)= 1; 
     end
end
A
% prosdiorismos tou dianysmatos b
% Symfwna me ton orismo dinw times sta dianysmata
for i=1:n-1
   b_y(i) =  (6/(h^2))* (y(i) - 2*y(i+1) + y(i+2)); % Dimiourgo ena dianysma by, n-1 orwn gia ton deikth TSE300
   b_z(i) =  (6/(h^2))* (z(i) - 2*z(i+1) + z(i+2)); % Dimiourgo ena dianysma bz, n-1 orwn gia ton deikth NASDAQ
end

%% (3) Ypologismos twn timwn mj, gia kathe synolo dedomenwn
% Sthn periptwsh ayth tha prepi na lisw to systhma Am=b

% Ypologizw ta m(2), m(3), ..., m(n) gia tou diktes TSE300 kai NASDAQ

tempy = A\b_y';  %deikth TSE300
tempz = A\b_z';  %deikth NASDAQ

% Ypologismos twn m(j), j=1,2,...,n,n+1. gia ton deikth TSE300
% Ypologismos twn m(j), j=1,2,...,n,n+1. gia ton deikth NASDAQ
for i=2:n
        my(i)=tempy(i-1);   %deikths TSE300
        mz(i)=tempz(i-1);   %deikths NASDAQ
end
my(1) = (2*my(2)) - my(3);  %deikths TSE300
mz(1) = (2*mz(2)) - mz(3);  %deikths NASDAQ
my(n+1) = (2*my(n)) - my(n-1);  %deikths TSE300
mz(n+1) = (2*mz(n)) - mz(n-1);  %deikths NASDAQ

my
mz

%% (4) ypologismos twn paramertwn ( a(j), b(j), c(j), d(j) ), twn
%% polywnymown kivikhs paremvolhs y = Sj(x)

ay = y(1:n);
az = z(1:n);
   
by = (y(2:n+1) - y(1:n))./h - h.*(my(2:n+1)+2*my(1:n))/6;
bz = (z(2:n+1) - z(1:n))./h - h.*(mz(2:n+1)+2*mz(1:n))/6;

   
cy = 0.5*my(1:n);
cz = 0.5*mz(1:n);

dy = (my(2:n+1) - my(1:n))./(6*h);
dz = (mz(2:n+1) - mz(1:n))./(6*h);

%% (5) Sxediasmos twn splines paremvolis sto diasthma[0,5] gia ta 2 synola
%% (6) kai Sxediasmos twn simiwn twn antoistixwn dedomenwn

xInt = 1 : 0.01 : 5;
for j = 1 : length(xInt)
	if xInt(j) ~= x(n+1)
        iInt(j) = sum(x <= xInt(j));
	else
		iInt(j) = n;
	end
end
yInt = ay(iInt) + by(iInt).*(xInt-x(iInt)) + cy(iInt).*(xInt - x(iInt)).^2 + dy(iInt).*(xInt - x(iInt)).^3;
plot(xInt,yInt,'g',x,y,'b*');   axis([0,5,7500,7950]);

pause;
zInt = az(iInt) + bz(iInt).*(xInt-x(iInt)) + cz(iInt).*(xInt - x(iInt)).^2 + dz(iInt).*(xInt - x(iInt)).^3;
plot(xInt,zInt,'g',x,z,'b*');   axis([0,5,1800,1870]);


% Oi times sta mesa ths tetarths
yInt = ay(iInt) + by(iInt).*(3.5-x(iInt)) + cy(iInt).*(3.5 - x(iInt)).^2 + dy(iInt).*(3.5 - x(iInt)).^3;
zInt = az(iInt) + bz(iInt).*(3.5-x(iInt)) + cz(iInt).*(3.5 - x(iInt)).^2 + dz(iInt).*(3.5 - x(iInt)).^3;

fprintf('I timh sta mesa ths tetarths:\n');
fprintf('Gia ton deikth TSE300 :');disp(sum(yInt)/length(yInt))
fprintf('\nGia ton deikth NASDAQ :');disp(sum(zInt)/length(zInt))
pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Ypologismos sta mesa ths tetarths me thn synarthsh spline tou Matlab
fprintf('Ypologimos sta mesa ths tetarths me thn synarthsh spline tou Matlab\n');
m_tetarths = spline(x,y,3.5);
disp(m_tetarths)
xx = 0:0.01:5;
yy = spline(x,y,xx);
plot(xx,yy,'g',x,y,'b*',3.5,m_tetarths,'r*');
      
    
